library(SDMTools)
species.data=read.csv('/home/jc112725/CARABIDS/CARABIDS/species_list_with_range1990.csv')
species=species.data$Species; species=as.vector(species)


#set working directory where the raster files will be pulled from
work.dir = '/home/jc112725/CARABIDS/models/'; setwd(work.dir)

#create empty object
richness= NULL
#make  aloop opening first asc file
for (spp in species){ cat(spp,'\n')
tasc=read.asc.gz(paste(work.dir,spp,"/summary/1990.realized.asc.gz", sep=""))
tasc[which(tasc>0)]=1 #set all suitable area to 1

if (length (richness)==0){richness=tasc} else {richness=richness+tasc}
}

#write.asc(richness,"/home/jc112725/CARABIDS/richness.asc")
max.richness = max(richness, na.rm=TRUE) #get the max richness value
pnts=cbind(x=c(146,146.15,146.15,146), y=c(-15.7,-15.7,-16.5,-16.5))
zlim=c(0,max.richness)
Colourmap= c("deepskyblue4",heat.colors(100)[seq(90,1)])

png ("/home/jc112725/CARABIDS/richness.png",bg= "white",height=800, width=400)
par(cex=1,mar=c(3,3,1,1))
         image(richness, col=Colourmap,zlim=zlim)
 legend.gradient(pnts,col=Colourmap,limits=c(0,max.richness), title='Species richness', cex=1)
   dev.off()






